The aim of this practical is to enhance your understanding of multiple imputation, in general. You will learn how to pool the results of analyses performed on multiply-imputed data, how to approach different types of data and how to avoid the pitfalls researchers may fall into. The main objective is to increase your knowledge and understanding on applications of multiple imputation.
Again, we start by loading (with library()) the necessary packages and fixing the random seed to allow for our outcomes to be replicable.
library(mice)
library(lme4)
library(mitml)
set.seed(123)
All the best,
popularity data setWe are going to work with the popularity data from Joop Hox (2010). The variables in this data set are described as follows:
| pupil | \(\quad\) Pupil number within class |
| class | \(\quad\) Class number |
| extrav | \(\quad\) Pupil extraversion |
| sex | \(\quad\) Pupil gender |
| texp | \(\quad\) Teacher experience (years) |
| popular | \(\quad\) Popularity as rated by the pupil’s peers |
| popteach | \(\quad\) Popularity as rated by the teacher |
A workspace with complete and incomplete versions of the popularity data can be obtained here or can be loaded into the Global Environment by running:
con <- url("https://www.gerkovink.com/mimp/popular.RData")
load(con)
This workspace contains several datasets and functions that, when loaded, are available to you in R. If you’d like to see what is in your Gloval Environment after importing the workspace: run the following code
ls()
## [1] "con" "icc"
## [3] "mice.impute.2lonly.mean2" "popMCAR"
## [5] "popMCAR2" "popNCR"
## [7] "popNCR2" "popNCR3"
## [9] "popular"
The popular data set is the complete version. The other popXXX sets are incomplete versions of this data set with varying levels of complexity. The icc() function is a generic function from package multilevel (function ICC1()).
Having the complete version of the data (the popular set) is a luxury that we do not have in practice. However, for educational purposes we can use this simulated set and its incomplete versions to inspect, evaluate and verify our imputation methodology and approach.
The dataset popNCR is a variation on the Hox (2010) data, where the missingness in the variables is either missing at random (MAR) or missing not at random (MNAR). We will be working with this popNCR set in this practical.
popNCR data1. Check with the functions head(), dim() - alternatively one could use nrow() and ncol() instead of dim() - and summary() how large the dataset is, of which variables the data frame consists and if there are missing values in a variable.
head(popNCR)
## pupil class extrav sex texp popular popteach
## 1 1 1 5 1 NA 6.3 NA
## 2 2 1 NA 0 24 4.9 NA
## 3 3 1 4 1 NA 5.3 6
## 4 4 1 3 <NA> NA 4.7 5
## 5 5 1 5 1 24 NA 6
## 6 6 1 NA 0 NA 4.7 5
dim(popNCR)
## [1] 2000 7
nrow(popNCR)
## [1] 2000
ncol(popNCR)
## [1] 7
summary(popNCR)
## pupil class extrav sex texp
## Min. : 1.00 17 : 26 Min. : 1.000 0 :661 Min. : 2.0
## 1st Qu.: 6.00 63 : 25 1st Qu.: 4.000 1 :843 1st Qu.: 7.0
## Median :11.00 10 : 24 Median : 5.000 NA's:496 Median :12.0
## Mean :10.65 15 : 24 Mean : 5.313 Mean :11.8
## 3rd Qu.:16.00 4 : 23 3rd Qu.: 6.000 3rd Qu.:16.0
## Max. :26.00 21 : 23 Max. :10.000 Max. :25.0
## (Other):1855 NA's :516 NA's :976
## popular popteach
## Min. :0.000 Min. : 1.000
## 1st Qu.:3.900 1st Qu.: 4.000
## Median :4.800 Median : 5.000
## Mean :4.829 Mean : 4.834
## 3rd Qu.:5.800 3rd Qu.: 6.000
## Max. :9.100 Max. :10.000
## NA's :510 NA's :528
The data set has 2000 rows and 7 columns (variables). The variables extrav, sex, texp, popular and popteach contain missings. About a quarter of these variables is missing, except for texp where 50 % is missing.
2. As we have seen before, the function md.pattern() is used to display all different missing data patterns. How many different missing data patterns are present in the popNCR dataframe and which patterns occur most frequently in the data?
md.pattern(popNCR, plot = FALSE)
## pupil class sex popular extrav popteach texp
## 308 1 1 1 1 1 1 1 0
## 279 1 1 1 1 1 1 0 1
## 110 1 1 1 1 1 0 1 1
## 115 1 1 1 1 1 0 0 2
## 114 1 1 1 1 0 1 1 1
## 98 1 1 1 1 0 1 0 2
## 33 1 1 1 1 0 0 1 2
## 24 1 1 1 1 0 0 0 3
## 119 1 1 1 0 1 1 1 1
## 113 1 1 1 0 1 1 0 2
## 50 1 1 1 0 1 0 1 2
## 75 1 1 1 0 1 0 0 3
## 29 1 1 1 0 0 1 1 2
## 21 1 1 1 0 0 1 0 3
## 2 1 1 1 0 0 0 1 3
## 14 1 1 1 0 0 0 0 4
## 102 1 1 0 1 1 1 1 1
## 89 1 1 0 1 1 1 0 2
## 25 1 1 0 1 1 0 1 2
## 29 1 1 0 1 1 0 0 3
## 85 1 1 0 1 0 1 1 2
## 56 1 1 0 1 0 1 0 3
## 9 1 1 0 1 0 0 1 3
## 14 1 1 0 1 0 0 0 4
## 19 1 1 0 0 1 1 1 2
## 27 1 1 0 0 1 1 0 3
## 13 1 1 0 0 1 0 1 3
## 11 1 1 0 0 1 0 0 4
## 4 1 1 0 0 0 1 1 3
## 9 1 1 0 0 0 1 0 4
## 2 1 1 0 0 0 0 1 4
## 2 1 1 0 0 0 0 0 5
## 0 0 496 510 516 528 976 3026
md.pattern(popNCR)
There are 32 unique patterns. The pattern where everything is observed and the pattern where only texp is missing occur most frequently.
If we omit texp, then half the patterns disappear:
md.pattern(popNCR[, -5], plot = FALSE)
## pupil class sex popular extrav popteach
## 587 1 1 1 1 1 1 0
## 225 1 1 1 1 1 0 1
## 212 1 1 1 1 0 1 1
## 57 1 1 1 1 0 0 2
## 232 1 1 1 0 1 1 1
## 125 1 1 1 0 1 0 2
## 50 1 1 1 0 0 1 2
## 16 1 1 1 0 0 0 3
## 191 1 1 0 1 1 1 1
## 54 1 1 0 1 1 0 2
## 141 1 1 0 1 0 1 2
## 23 1 1 0 1 0 0 3
## 46 1 1 0 0 1 1 2
## 24 1 1 0 0 1 0 3
## 13 1 1 0 0 0 1 3
## 4 1 1 0 0 0 0 4
## 0 0 496 510 516 528 2050
md.pattern(popNCR[, -5])
Without the fifth column [, -5] which is texp, there are only 16 patterns.
popular3. Let’s focus more precisely on the missing data patterns. Does the missing data of popular depend on popteach? One could for example check this by making a histogram of popteach separately for the pupils with known popularity and missing popularity.
In R the missingness indicator
is.na(popNCR$popular)
is a dummy variable of the same length as popular with value 0 (FALSE) for observed pupil popularity and 1 (TRUE) for missing pupil popularity. A histogram can be made with the function histogram(). The code for a conditional histogram of popteach given the missingness indicator for popular is
histogram(~ popteach | is.na(popular), data=popNCR)
We do see that the histogram for the missing popular (TRUE) is further to the right than the histogram for observed popular (FALSE). This would indicate a right-tailed MAR missingness. In fact this is exactly what happens, because I (Gerko) created the missingness in these data. However, we can make it observable by examining the relations between the missingness in popular and the observed data in popteach. This is a clear indication that the missingness in popular depends on the observed values for popteach. This finding makes the MCAR assumption improbable.
4. Does the missingness of the other incomplete variables depend on popteach? If yes, what is the direction of the relation?
histogram(~ popteach | is.na(sex), data = popNCR)
There seems to be a left-tailed relation between popteach and the missingness in sex.
histogram(~ popteach | is.na(extrav), data = popNCR)
There also seems to be a left-tailed relation between popteach and the missingness in extrav.
histogram(~ popteach | is.na(texp), data = popNCR)
There seems to be no observable relation between popteach and the missingness in texp. It might be random or even MNAR.
5. Find out if the missingness in teacher popularity depends on pupil popularity.
histogram(~ popular | is.na(popteach), data = popNCR)
Yes: there is a dependency. The relation seems to be right-tailed.
We now know that some variables (columns) seem to be related to the missingness in other variables. This is a clear indication that for the relevant variables, the predictor relations during multiple imputation should be set.
6. Have a look at the intraclass correlation (ICC) for the incomplete variables popular, popteach and texp.
icc(aov(popular ~ as.factor(class), data = popNCR))
## [1] 0.328007
icc(aov(popteach ~ class, data = popNCR))
## [1] 0.3138658
icc(aov(texp ~ class, data = popNCR))
## [1] 1
Please note that the function icc() comes from the package multilevel (function ICC1()), but is included in the workspace popular.RData. Make a note of the ICC’s, you’ll need them later in Practical 3.
7. Do you think it is necessary to take the multilevel structure into account?
YES! There is a strong cluster structure going on. If we ignore the clustering in our imputation model, we may run into invalid inference. To stay as close to the true data model, we must take the cluster structure into account during imputation.
popNCR data set8. Impute the popNCR dataset with mice using imputation method norm for popular, popteach, texp and extrav. Exclude class as a predictor for all variables. Name the multiply imputed data set (mids) resulting from the mice call imp1. Use \(m=10\) imputations and \(maxit = 15\) iterations.
#create adjusted imputation method vector
meth <- make.method(popNCR)
meth[c("popular", "popteach", "texp", "extrav")] <- "norm"
meth
## pupil class extrav sex texp popular popteach
## "" "" "norm" "logreg" "norm" "norm" "norm"
#create adjusted predictorMatrix
pred <- make.predictorMatrix(popNCR)
pred[, "class"] <- 0
pred[, "pupil"] <- 0
pred
## pupil class extrav sex texp popular popteach
## pupil 0 0 1 1 1 1 1
## class 0 0 1 1 1 1 1
## extrav 0 0 0 1 1 1 1
## sex 0 0 1 0 1 1 1
## texp 0 0 1 1 0 1 1
## popular 0 0 1 1 1 0 1
## popteach 0 0 1 1 1 1 0
#generate imputations
imp1 <- mice(popNCR, m = 10, maxit = 15,
method = meth,
predictorMatrix = pred,
print = FALSE)
9. Inspect convergence of the mice algorithm by studying the traceplots
plot(imp1)
Convergence has been reached for some variables, but is not convincing for all. However, we’ve not done any justice to the structure of the data - we’ve imputed the data as a flat file without paying any attention to the clustering in the data. The imputation model can be greatly improved, which we will see in the next practical.
9. Run the following model on the imputed data, the incomplete data and compare the model to the ‘true’ data
library(lme4)
fit.true <- with(popular, lmer(popular ~ 1 + (1|class)))
summary(fit.true)
## Linear mixed model fit by REML ['lmerMod']
## Formula: popular ~ 1 + (1 | class)
##
## REML criterion at convergence: 6330.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5655 -0.6975 0.0020 0.6758 3.3175
##
## Random effects:
## Groups Name Variance Std.Dev.
## class (Intercept) 0.7021 0.8379
## Residual 1.2218 1.1053
## Number of obs: 2000, groups: class, 100
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.07786 0.08739 58.1
fit.true <- with(popNCR, lmer(popular ~ 1 + (1|class)))
summary(fit.true)
## Linear mixed model fit by REML ['lmerMod']
## Formula: popular ~ 1 + (1 | class)
##
## REML criterion at convergence: 4727.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2577 -0.6985 -0.0257 0.6649 3.2828
##
## Random effects:
## Groups Name Variance Std.Dev.
## class (Intercept) 0.6298 0.7936
## Residual 1.2078 1.0990
## Number of obs: 1490, groups: class, 100
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.91988 0.08452 58.21
fit <- with(imp1, lmer(popular ~ 1 + (1|class)))
pool <- pool(fit)
summary(pool)
## estimate std.error statistic df p.value
## (Intercept) 5.013788 0.0772285 64.92147 1731.811 0
We may obtain the variance components by the testEstimates() function from package mitml:
library(mitml)
testEstimates(as.mitml.result(fit), var.comp = TRUE)$var.comp
## Estimate
## Intercept~~Intercept|class 0.5153594
## Residual~~Residual 1.3231304
## ICC|class 0.2803254
There are quite some differences between the models. With multiple imputation, we are able to correct the inference quite to some extend (e.g. the fixed effects), but we miss out on correcting the inference on the level of the random effects. Then again, we completely ignored the clustering structure in our data.
10. Compare the means of the variables in the first imputed dataset and in the incomplete dataset.
summary(complete(imp1))
## pupil class extrav sex
## Min. : 1.00 17 : 26 Min. : 0.3408 0: 973
## 1st Qu.: 6.00 63 : 25 1st Qu.: 4.0000 1:1027
## Median :11.00 10 : 24 Median : 5.0000
## Mean :10.65 15 : 24 Mean : 5.2210
## 3rd Qu.:16.00 4 : 23 3rd Qu.: 6.0000
## Max. :26.00 21 : 23 Max. :10.0000
## (Other):1855
## texp popular popteach
## Min. :-7.887 Min. :0.000 Min. : 1.000
## 1st Qu.: 8.000 1st Qu.:4.044 1st Qu.: 4.000
## Median :12.976 Median :5.000 Median : 5.000
## Mean :12.612 Mean :4.986 Mean : 5.011
## 3rd Qu.:17.000 3rd Qu.:5.900 3rd Qu.: 6.000
## Max. :33.472 Max. :9.100 Max. :10.000
##
summary(popNCR)
## pupil class extrav sex texp
## Min. : 1.00 17 : 26 Min. : 1.000 0 :661 Min. : 2.0
## 1st Qu.: 6.00 63 : 25 1st Qu.: 4.000 1 :843 1st Qu.: 7.0
## Median :11.00 10 : 24 Median : 5.000 NA's:496 Median :12.0
## Mean :10.65 15 : 24 Mean : 5.313 Mean :11.8
## 3rd Qu.:16.00 4 : 23 3rd Qu.: 6.000 3rd Qu.:16.0
## Max. :26.00 21 : 23 Max. :10.000 Max. :25.0
## (Other):1855 NA's :516 NA's :976
## popular popteach
## Min. :0.000 Min. : 1.000
## 1st Qu.:3.900 1st Qu.: 4.000
## Median :4.800 Median : 5.000
## Mean :4.829 Mean : 4.834
## 3rd Qu.:5.800 3rd Qu.: 6.000
## Max. :9.100 Max. :10.000
## NA's :510 NA's :528
11. The missingness in texp is MNAR: higher values for texp have a larger probability to be missing. Can you see this in the imputed data? Do you think this is a problem?
Yes, we can see this in the imputed data: teacher experience increases slightly after imputation. However, texp is the same for all pupils in a class. But not all pupils have this information recorded (as if some pupils did not remember, or were not present during data collection). This is not a problem, because as long as at least one pupil in each class has teacher experience recorded, we can deductively impute the correct (i.e. true) value for every pupil in the class.
12. Compare the ICC’s of the variables in the first imputed dataset with those in the incomplete dataset (use popular, popteach and texp). Make a notation of the ICC’s after imputation.
data.frame(vars = names(popNCR[c(6, 7, 5)]),
observed = c(icc(aov(popular ~ class, popNCR)),
icc(aov(popteach ~ class, popNCR)),
icc(aov(texp ~ class, popNCR))),
norm = c(icc(aov(popular ~ class, complete(imp1))),
icc(aov(popteach ~ class, complete(imp1))),
icc(aov(texp ~ class, complete(imp1)))))
## vars observed norm
## 1 popular 0.3280070 0.2785369
## 2 popteach 0.3138658 0.2542062
## 3 texp 1.0000000 0.4324512
13. Now impute the popNCR dataset again with mice using imputation method norm for popular, popteach, texp and extrav, but now include class as a predictor for all variables. Call the mids-object imp2.
pred <- make.predictorMatrix(popNCR)
pred[, "pupil"] <- 0
imp2 <- mice(popNCR, meth = meth, pred = pred, print = FALSE)
## Warning: Number of logged events: 90
We exclude pupil here to avoid overfitting our data with the pupil identifier. Including pupil would result in a model with zero residual variance.
14. Compare the ICC’s of the variables in the first imputed dataset from imp2 with those of imp1 and the incomplete dataset (use popular, popteach and texp). Make a notation of the ICC’s after imputation.
data.frame(vars = names(popNCR[c(6, 7, 5)]),
observed = c(icc(aov(popular ~ class, popNCR)),
icc(aov(popteach ~ class, popNCR)),
icc(aov(texp ~ class, popNCR))),
norm = c(icc(aov(popular ~ class, complete(imp1))),
icc(aov(popteach ~ class, complete(imp1))),
icc(aov(texp ~ class, complete(imp1)))),
normclass = c(icc(aov(popular ~ class, complete(imp2))),
icc(aov(popteach ~ class, complete(imp2))),
icc(aov(texp ~ class, complete(imp2)))))
## vars observed norm normclass
## 1 popular 0.3280070 0.2785369 0.3670845
## 2 popteach 0.3138658 0.2542062 0.3347250
## 3 texp 1.0000000 0.4324512 1.0000000
By simply forcing the algorithm to use the class variable during estimation we adopt a fixed effects approach. This conforms to formulating seperate regression models for each class and imputing within classes from these models.
15. Inspect the trace lines for the variables popular, texp and extrav.
plot(imp2, c("popular", "texp", "popteach"))
It seems not quite OK. For example, there are clear trends for popteach and some of the streams do not even intermingle between iterations 4 and 5 for texp. Adding another 20 iterations gives a more convincing graph:
imp2b <- mice.mids(imp2, maxit = 20, print = FALSE)
plot(imp2b, c("popular", "texp", "popteach"))
We can use function mice.mids() to continue from the state where the object imp2 stopped. This is fortunate, as we do not have to re-run the first maxit = 5 default iterations.
16. Plot the densities of the observed and imputed data (use imp2) with the function densityplot().
To obtain all densities of the different imputed datasets use
densityplot(imp2)
It is always wise to ask yourself the following questions when attacking a multilevel missing data problem - see Table 7.1 in Van Buuren (2018):
Q1. Will the complete-data model include random slopes?: Thus far we have not formulated any complicated analysis model, but if our analysis model would include random slopes than we should definitely account for those terms in our imputation model. For now, let’s assume that our analysis model remains: lme4::lmer(popular ~ (1 | class))
Q2. Will the data contain systematically missing values? With systematically missing data there are no observed values in the cluster
We can simply verify this by asking the length of the complete cases for each cluster
aggregate(. ~ class, data = popNCR, FUN=function(x) length(is.na(x)))
## class pupil extrav sex texp popular popteach
## 1 1 4 4 4 4 4 4
## 2 2 4 4 4 4 4 4
## 3 3 4 4 4 4 4 4
## 4 4 1 1 1 1 1 1
## 5 5 7 7 7 7 7 7
## 6 6 5 5 5 5 5 5
## 7 8 5 5 5 5 5 5
## 8 9 5 5 5 5 5 5
## 9 11 2 2 2 2 2 2
## 10 12 3 3 3 3 3 3
## 11 13 2 2 2 2 2 2
## 12 14 3 3 3 3 3 3
## 13 15 4 4 4 4 4 4
## 14 16 3 3 3 3 3 3
## 15 17 3 3 3 3 3 3
## 16 18 3 3 3 3 3 3
## 17 19 4 4 4 4 4 4
## 18 20 2 2 2 2 2 2
## 19 21 3 3 3 3 3 3
## 20 22 4 4 4 4 4 4
## 21 23 2 2 2 2 2 2
## 22 24 6 6 6 6 6 6
## 23 25 1 1 1 1 1 1
## 24 26 3 3 3 3 3 3
## 25 27 2 2 2 2 2 2
## 26 28 5 5 5 5 5 5
## 27 30 3 3 3 3 3 3
## 28 31 1 1 1 1 1 1
## 29 32 7 7 7 7 7 7
## 30 33 2 2 2 2 2 2
## 31 34 6 6 6 6 6 6
## 32 35 3 3 3 3 3 3
## 33 36 4 4 4 4 4 4
## 34 37 1 1 1 1 1 1
## 35 38 4 4 4 4 4 4
## 36 39 2 2 2 2 2 2
## 37 40 3 3 3 3 3 3
## 38 41 1 1 1 1 1 1
## 39 42 6 6 6 6 6 6
## 40 43 5 5 5 5 5 5
## 41 44 3 3 3 3 3 3
## 42 45 8 8 8 8 8 8
## 43 46 2 2 2 2 2 2
## 44 47 1 1 1 1 1 1
## 45 48 3 3 3 3 3 3
## 46 49 4 4 4 4 4 4
## 47 50 1 1 1 1 1 1
## 48 51 2 2 2 2 2 2
## 49 52 7 7 7 7 7 7
## 50 53 1 1 1 1 1 1
## 51 54 4 4 4 4 4 4
## 52 55 4 4 4 4 4 4
## 53 56 6 6 6 6 6 6
## 54 57 3 3 3 3 3 3
## 55 58 2 2 2 2 2 2
## 56 59 7 7 7 7 7 7
## 57 60 7 7 7 7 7 7
## 58 61 1 1 1 1 1 1
## 59 62 2 2 2 2 2 2
## 60 63 5 5 5 5 5 5
## 61 65 1 1 1 1 1 1
## 62 66 1 1 1 1 1 1
## 63 67 1 1 1 1 1 1
## 64 68 3 3 3 3 3 3
## 65 69 2 2 2 2 2 2
## 66 70 1 1 1 1 1 1
## 67 71 2 2 2 2 2 2
## 68 72 1 1 1 1 1 1
## 69 73 8 8 8 8 8 8
## 70 74 2 2 2 2 2 2
## 71 75 1 1 1 1 1 1
## 72 76 2 2 2 2 2 2
## 73 77 5 5 5 5 5 5
## 74 78 2 2 2 2 2 2
## 75 79 1 1 1 1 1 1
## 76 80 2 2 2 2 2 2
## 77 81 6 6 6 6 6 6
## 78 82 6 6 6 6 6 6
## 79 84 1 1 1 1 1 1
## 80 85 1 1 1 1 1 1
## 81 86 3 3 3 3 3 3
## 82 87 9 9 9 9 9 9
## 83 88 4 4 4 4 4 4
## 84 89 1 1 1 1 1 1
## 85 90 6 6 6 6 6 6
## 86 91 3 3 3 3 3 3
## 87 92 3 3 3 3 3 3
## 88 93 3 3 3 3 3 3
## 89 94 2 2 2 2 2 2
## 90 95 3 3 3 3 3 3
## 91 96 4 4 4 4 4 4
## 92 97 2 2 2 2 2 2
## 93 99 1 1 1 1 1 1
## 94 100 4 4 4 4 4 4
None of the clusters class are completely unobserved.
Q3. Will the distribution of the residuals be non-normal?
fit <- with(popNCR, lmer(popular ~ (1 | class)))
res <- residuals(fit)
hist(res, main = "Residuals")
The residuals seem to be quite normally distributed.
Q4. Will the error variance differ over clusters?
dotplot(ranef(fit, condVar = TRUE))
## $class
The error variance seems quite constant from a quick visual inspection
Q5. Will there be small clusters?
aggregate(pupil ~ class, data = popNCR, length)
## class pupil
## 1 1 20
## 2 2 20
## 3 3 18
## 4 4 23
## 5 5 21
## 6 6 20
## 7 7 21
## 8 8 20
## 9 9 20
## 10 10 24
## 11 11 22
## 12 12 17
## 13 13 20
## 14 14 17
## 15 15 24
## 16 16 17
## 17 17 26
## 18 18 21
## 19 19 20
## 20 20 20
## 21 21 23
## 22 22 20
## 23 23 22
## 24 24 19
## 25 25 20
## 26 26 21
## 27 27 20
## 28 28 20
## 29 29 17
## 30 30 17
## 31 31 19
## 32 32 23
## 33 33 19
## 34 34 22
## 35 35 20
## 36 36 19
## 37 37 19
## 38 38 17
## 39 39 16
## 40 40 20
## 41 41 19
## 42 42 21
## 43 43 17
## 44 44 19
## 45 45 19
## 46 46 23
## 47 47 17
## 48 48 19
## 49 49 20
## 50 50 20
## 51 51 23
## 52 52 23
## 53 53 19
## 54 54 19
## 55 55 21
## 56 56 19
## 57 57 18
## 58 58 20
## 59 59 17
## 60 60 21
## 61 61 19
## 62 62 20
## 63 63 25
## 64 64 19
## 65 65 22
## 66 66 20
## 67 67 21
## 68 68 17
## 69 69 18
## 70 70 21
## 71 71 22
## 72 72 21
## 73 73 20
## 74 74 17
## 75 75 17
## 76 76 18
## 77 77 19
## 78 78 18
## 79 79 20
## 80 80 20
## 81 81 21
## 82 82 21
## 83 83 18
## 84 84 21
## 85 85 22
## 86 86 22
## 87 87 21
## 88 88 23
## 89 89 18
## 90 90 23
## 91 91 18
## 92 92 17
## 93 93 20
## 94 94 16
## 95 95 21
## 96 96 21
## 97 97 21
## 98 98 21
## 99 99 23
## 100 100 20
We can see that the average cluster size is approximately 20 pupils in every cluster with no cluster having fewer than 16 pupils.
Q6. Will there be a small number of clusters? NO, there are 100 clusters
Q7. Will the complete-data model have cross-level interactions? No
Q8. Will the dataset be very large? No, only 2000 cases over 100 clusters
** End of intermezzo**
norm with pmm17. Impute the popNCR data once more where you use predictive mean matching and exclude only pupil as a predictor. Name the object imp4.
imp4 <- mice(popNCR[, -1])
##
## iter imp variable
## 1 1 extrav sex texp popular popteach
## 1 2 extrav sex texp popular popteach
## 1 3 extrav sex texp popular popteach
## 1 4 extrav sex texp popular popteach
## 1 5 extrav sex texp popular popteach
## 2 1 extrav sex texp popular popteach
## 2 2 extrav sex texp popular popteach
## 2 3 extrav sex texp popular popteach
## 2 4 extrav sex texp popular popteach
## 2 5 extrav sex texp popular popteach
## 3 1 extrav sex texp popular popteach
## 3 2 extrav sex texp popular popteach
## 3 3 extrav sex texp popular popteach
## 3 4 extrav sex texp popular popteach
## 3 5 extrav sex texp popular popteach
## 4 1 extrav sex texp popular popteach
## 4 2 extrav sex texp popular popteach
## 4 3 extrav sex texp popular popteach
## 4 4 extrav sex texp popular popteach
## 4 5 extrav sex texp popular popteach
## 5 1 extrav sex texp popular popteach
## 5 2 extrav sex texp popular popteach
## 5 3 extrav sex texp popular popteach
## 5 4 extrav sex texp popular popteach
## 5 5 extrav sex texp popular popteach
## Warning: Number of logged events: 90
18. Plot again the densities of the observed and imputed data with the function densityplot(), but now use imp4. Is there a difference between the imputations obtained with pmm and norm and can you explain this?
densityplot(imp4)
Yes, pmm samples from the observed values and this clearly shows: imputations follow the shape of the observed data.
19. Compare the ICC’s of the variables in the first imputed dataset from imp4 with those of imp1, imp2 and the incomplete dataset (use popular, popteach and texp).
See Exercise 12 for the solution.
20. Finally, compare the ICC’s of the imputations to the ICC’s in the original data. The original data can be found in dataset popular. What do you conclude?
data.frame(vars = names(popNCR[c(6, 7, 5)]),
observed = c(icc(aov(popular ~ class, popNCR)),
icc(aov(popteach ~ class, popNCR)),
icc(aov(texp ~ class, popNCR))),
norm = c(icc(aov(popular ~ class, complete(imp1))),
icc(aov(popteach ~ class, complete(imp1))),
icc(aov(texp ~ class, complete(imp1)))),
normclass = c(icc(aov(popular ~ class, complete(imp2))),
icc(aov(popteach ~ class, complete(imp2))),
icc(aov(texp ~ class, complete(imp2)))),
pmm = c(icc(aov(popular ~ class, complete(imp4))),
icc(aov(popteach ~ class, complete(imp4))),
icc(aov(texp ~ class, complete(imp4)))),
orig = c(icc(aov(popular ~ as.factor(class), popular)),
icc(aov(popteach ~ as.factor(class), popular)),
icc(aov(texp ~ as.factor(class), popular))))
## vars observed norm normclass pmm orig
## 1 popular 0.3280070 0.2785369 0.3670845 0.3701673 0.3629933
## 2 popteach 0.3138658 0.2542062 0.3347250 0.3433643 0.3414766
## 3 texp 1.0000000 0.4324512 1.0000000 1.0000000 1.0000000
Note: these display only the first imputed data set.
Mice includes several imputation methods for imputing multilevel data:
lme4::lmer()lme4::glmer()The latter two methods aggregate level 1 variables at level 2, but in combination with mice.impute.2l.pan, allow switching regression imputation between level 1 and level 2 as described in Yucel (2008) or Gelman and Hill (2006, p. 541). For more information on these imputation methods see the help.
There are also many useful imputation routines available from packages miceadds and micemd. See Tables 7.2, 7.3 and 7.4 in Van Buuren (2018) for a detailed overview of these functions.
21. Impute the variable popular by means of 2l.norm. Use dataset popNCR2.
ini <- mice(popNCR2, maxit = 0)
pred <- ini$pred
pred["popular", ] <- c(0, -2, 2, 2, 2, 0, 2)
In the predictor matrix, -2 denotes the class variable, a value 1 indicates a fixed effect and a value 2 indicates a random effect. However, the currently implemented algorithm does not handle predictors that are specified as fixed effects (type = 1). When using mice.impute.2l.norm(), the current advice is to specify all predictors as random effects (type = 2).
meth <- ini$meth
meth <- c("", "", "", "", "", "2l.norm", "")
imp5 <- mice(popNCR2, pred = pred, meth=meth, print = FALSE)
22. Inspect the imputations. Did the algorithm converge?
densityplot(imp5, ~popular, ylim = c(0, 0.35), xlim = c(-1.5, 10))
densityplot(imp4, ~popular, ylim = c(0, 0.35), xlim = c(-1.5, 10))
The imputations generated with 2l.norm are very similar to the ones obtained by pmm with class as a fixed effect. If we plot the first imputed dataset from imp4 and imp5 against the original (true) data:
plot(density(popular$popular)) #true data
lines(density(complete(imp5)$popular), col = "red", lwd = 2) #2l.norm
lines(density(complete(imp4)$popular), col = "green", lwd = 2) #PMM
We can see that the imputations are very similar. When studying the convergence
plot(imp5)
we conclude that it may be wise to run additional iterations. Convergence is not apparent from this plot.
imp5.b <- mice.mids(imp5, maxit = 10, print = FALSE)
plot(imp5.b)
After running another 10 iterations, convergence is more convincing.
23. In the original data, the group variances for popular are homogeneous. Use 2l.pan to impute the variable popular in dataset popNCR2. Inspect the imputations. Did the algorithm converge?
ini <- mice(popNCR2, maxit = 0)
pred <- ini$pred
pred["popular", ] <- c(0, -2, 2, 2, 1, 0, 2)
meth <- ini$meth
meth <- c("", "", "", "", "", "2l.pan", "")
imp6 <- mice(popNCR2, pred = pred, meth = meth, print = FALSE)
Let us create the densityplot for imp6
densityplot(imp6, ~popular, ylim = c(0, 0.35), xlim = c(-1.5, 10))
and compare it to the one for imp4
densityplot(imp4, ~popular, ylim = c(0, 0.35), xlim = c(-1.5, 10))
If we plot the first imputed dataset from both objects against the original (true) density, we obtain the following plot:
plot(density(popular$popular), main = "black = truth | green = PMM | red = 2l.pan") #
lines(density(complete(imp6)$popular), col = "red", lwd = 2) #2l.pan
lines(density(complete(imp4)$popular), col = "green", lwd = 2) #PMM
We can see that the imputations are very similar. When studying the convergence
plot(imp6)
we conclude that it may be wise to run additional iterations. Convergence is not apparent from this plot.
imp6.b <- mice.mids(imp5, maxit = 10, print = FALSE)
plot(imp6.b)
Again, after running another 10 iterations, convergence is more convincing.
popNCR324. Now inspect dataset popNCR3 and impute the incomplete variables according to the following imputation methods:
| Variable | Method |
|---|---|
| extrav | 2l.norm |
| texp | 2lonly.mean |
| sex | 2l.bin |
| popular | 2l.pan |
| popteach | 2l.lmer |
ini <- mice(popNCR3, maxit = 0)
pred <- ini$pred
pred["extrav", ] <- c(0, -2, 0, 2, 2, 2, 2) #2l.norm
pred["sex", ] <- c(0, -2, 1, 0, 1, 1, 1) #2l.bin
pred["texp", ] <- c(0, -2, 1, 1, 0, 1, 1) #2lonly.mean
pred["popular", ] <- c(0, -2, 2, 2, 1, 0, 2) #2l.pan
pred["popteach", ] <- c(0, -2, 2, 2, 1, 2, 0) #2l.lmer
meth <- ini$meth
meth <- c("", "", "2l.norm", "2l.bin", "2lonly.mean", "2l.pan", "2l.lmer")
imp7 <- mice(popNCR3, pred = pred, meth = meth, print = FALSE)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
fit <- with(imp7, lmer(popular ~ 1 + (1|class)))
testEstimates(as.mitml.result(fit), var.comp = TRUE)
##
## Call:
##
## testEstimates(model = as.mitml.result(fit), var.comp = TRUE)
##
## Final parameter estimates and inferences obtained from 5 imputed data sets.
##
## Estimate Std.Error t.value df P(>|t|) RIV FMI
## (Intercept) 5.043 0.088 57.591 53859.428 0.000 0.009 0.009
##
## Estimate
## Intercept~~Intercept|class 0.699
## Residual~~Residual 1.218
## ICC|class 0.364
##
## Unadjusted hypothesis test as appropriate in larger samples.
25. Evaluate the imputations by means of convergence, distributions and plausibility.
densityplot(imp7)
stripplot(imp7)
Given what we know about the missingness, the imputed densities look very reasonable.
plot(imp7)
Convergence has not yet been reached. more iterations are advisable.
pmm imputation with class as a predictor26 . Repeat the same imputations as in the previous step, but now use pmm for all non-binary incomplete variables. Include class as a predictor.
pmmdata <- popNCR3
pmmdata$class <- as.factor(popNCR3$class)
imp8 <- mice(pmmdata, m = 5, print = FALSE)
## Warning: Number of logged events: 90
A warning is printed that there are 90 logged events. When we inspect
imp8$loggedEvents
## it im dep meth out
## 1 1 1 popular pmm texp
## 2 1 1 popteach pmm texp
## 3 1 2 popular pmm texp
## 4 1 2 popteach pmm texp
## 5 1 3 popular pmm texp
## 6 1 3 popteach pmm texp
## 7 1 4 popular pmm texp
## 8 1 4 popteach pmm texp
## 9 1 5 popular pmm texp
## 10 1 5 popteach pmm texp
## 11 2 1 extrav pmm texp
## 12 2 1 sex logreg texp
## 13 2 1 popular pmm texp
## 14 2 1 popteach pmm texp
## 15 2 2 extrav pmm texp
## 16 2 2 sex logreg texp
## 17 2 2 popular pmm texp
## 18 2 2 popteach pmm texp
## 19 2 3 extrav pmm texp
## 20 2 3 sex logreg texp
## 21 2 3 popular pmm texp
## 22 2 3 popteach pmm texp
## 23 2 4 extrav pmm texp
## 24 2 4 sex logreg texp
## 25 2 4 popular pmm texp
## 26 2 4 popteach pmm texp
## 27 2 5 extrav pmm texp
## 28 2 5 sex logreg texp
## 29 2 5 popular pmm texp
## 30 2 5 popteach pmm texp
## 31 3 1 extrav pmm texp
## 32 3 1 sex logreg texp
## 33 3 1 popular pmm texp
## 34 3 1 popteach pmm texp
## 35 3 2 extrav pmm texp
## 36 3 2 sex logreg texp
## 37 3 2 popular pmm texp
## 38 3 2 popteach pmm texp
## 39 3 3 extrav pmm texp
## 40 3 3 sex logreg texp
## 41 3 3 popular pmm texp
## 42 3 3 popteach pmm texp
## 43 3 4 extrav pmm texp
## 44 3 4 sex logreg texp
## 45 3 4 popular pmm texp
## 46 3 4 popteach pmm texp
## 47 3 5 extrav pmm texp
## 48 3 5 sex logreg texp
## 49 3 5 popular pmm texp
## 50 3 5 popteach pmm texp
## 51 4 1 extrav pmm texp
## 52 4 1 sex logreg texp
## 53 4 1 popular pmm texp
## 54 4 1 popteach pmm texp
## 55 4 2 extrav pmm texp
## 56 4 2 sex logreg texp
## 57 4 2 popular pmm texp
## 58 4 2 popteach pmm texp
## 59 4 3 extrav pmm texp
## 60 4 3 sex logreg texp
## 61 4 3 popular pmm texp
## 62 4 3 popteach pmm texp
## 63 4 4 extrav pmm texp
## 64 4 4 sex logreg texp
## 65 4 4 popular pmm texp
## 66 4 4 popteach pmm texp
## 67 4 5 extrav pmm texp
## 68 4 5 sex logreg texp
## 69 4 5 popular pmm texp
## 70 4 5 popteach pmm texp
## 71 5 1 extrav pmm texp
## 72 5 1 sex logreg texp
## 73 5 1 popular pmm texp
## 74 5 1 popteach pmm texp
## 75 5 2 extrav pmm texp
## 76 5 2 sex logreg texp
## 77 5 2 popular pmm texp
## 78 5 2 popteach pmm texp
## 79 5 3 extrav pmm texp
## 80 5 3 sex logreg texp
## 81 5 3 popular pmm texp
## 82 5 3 popteach pmm texp
## 83 5 4 extrav pmm texp
## 84 5 4 sex logreg texp
## 85 5 4 popular pmm texp
## 86 5 4 popteach pmm texp
## 87 5 5 extrav pmm texp
## 88 5 5 sex logreg texp
## 89 5 5 popular pmm texp
## 90 5 5 popteach pmm texp
We find that texp has been excluded as a predictor in 90 instances. This is to be expected because when class is added as a factor (categorical variable) to the model, a seperate model will be fitted for each class. In each of these models, observed texp is a constant and, hence, will automatically be removed by mice to avoid estimation problems because of a redundant parameter.
With pmm, the imputations are very similar and conform to the shape of the observed data.
densityplot(imp8)
stripplot(imp8)
When looking at the convergence of pmm, more iterations are advisable:
plot(imp8)
fit <- with(imp8, lmer(popular ~ 1 + (1|class)))
testEstimates(as.mitml.result(fit), var.comp = TRUE)
##
## Call:
##
## testEstimates(model = as.mitml.result(fit), var.comp = TRUE)
##
## Final parameter estimates and inferences obtained from 5 imputed data sets.
##
## Estimate Std.Error t.value df P(>|t|) RIV FMI
## (Intercept) 5.047 0.088 57.136 21545.583 0.000 0.014 0.014
##
## Estimate
## Intercept~~Intercept|class 0.710
## Residual~~Residual 1.188
## ICC|class 0.374
##
## Unadjusted hypothesis test as appropriate in larger samples.
If the primary interest is on the fixed effects, using pmm in a fixed effects setup where the class is added as a cluster dummy variable, may be an easily implementable alternative to running a full multilevel model: see e.g. Vink, Lazendic & Van Buuren, 2015. In our example, when compared to the true data set, parameters for the model popular ~ 1 + (1|class) are unbiased. However, some authors have reported that the bias in random slopes and variance components can be substantial. If the focus lies on the random components in the model, a parametric full multilevel imputation would be more efficient if it is properly fitted. See Drechsler (2015), Lüdtke, Robitzsch, and Grund (2017) and Speidel, Drechsler, and Sakshaug (2017) for more detail on using cluster dummies in multilevel imputation and estimation.
There are ways to ensure that imputations are not just ‘guesses of unobserved values’. Imputations can be checked by using a standard of reasonability - see Aboyami, Gelman, and Levy (2010) for a wonderful introduction to diagnostics for multiple imputations. We are able to check the differences between observed and imputed values, the differences between their distributions as well as the distribution of the completed data as a whole. If we do this, we can see whether imputations make sense in the context of the problem being studied.
Abayomi, K. , Gelman, A. and Levy, M. (2008), Diagnostics for multivariate imputations. Journal of the Royal Statistical Society: Series C (Applied Statistics), 57: 273-291. Article
Drechsler, J. (2015). Multiple Imputation of Multilevel Missing Data: Rigor Versus Simplicity. Journal of Educational and Behavioral Statistics 40 (1): 69–95. Article
Gelman, A., & Hill, J. (2006). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.
Hox, J. J., Moerbeek, M., & van de Schoot, R. (2010). Multilevel analysis: Techniques and applications. Routledge.
Lüdtke, O., Robitzsch, A., & Grund, S. (2017). Multiple imputation of missing data in multilevel designs: A comparison of different strategies. Psychological Methods, 22(1), 141-165. Article
Speidel, M., J. Drechsler, and J. W. Sakshaug. (2017). Biases in Multilevel Analyses Caused by Cluster-Specific Fixed-Effects Imputation. Behavior Research Methods, 1–17. Article
Van Buuren, S. Flexible imputation of missing data. Chapman & Hall/CRC. 2018. online book
Vink, G., Lazendic, G., Van Buuren, S. (2015). Psychological Test and Assessment Modeling, volume 57, issue 4, pp. 577 - 594. Article
Yucel, R. M. (2008). Multiple imputation inference for multivariate multilevel continuous data with ignorable non-response. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 366(1874), 2389-2403. Article
- End of practical